#include "mcmcPLE.h"

////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
// for saving histogram or responses for one value at a certain position
MCMC_RespChange::MCMC_RespChange()
{
	num_resp = 0;
	pdata = NULL;
}

MCMC_RespChange::~MCMC_RespChange()
{
	num_resp = 0;
	if (pdata != NULL)
		delete []pdata;
	pdata = NULL;
}

void MCMC_RespChange::AllocateSpace(const int size)
{
	if (size!=num_resp)
	{
		if (pdata != NULL)
			delete []pdata;
		pdata = NULL;

		num_resp = size;
		if (num_resp>0)
			pdata = new char[size];
	}
}

void MCMC_RespChange::SaveToFile(ofstream &file)
{
	file.write((char*)&num_resp, sizeof(num_resp));
	if (num_resp>0)
		file.write((char*)pdata, num_resp*sizeof(char));
}

void MCMC_RespChange::ReadFromFile(ifstream &file)
{
	int num_temp;

	file.read((char*)&num_temp, sizeof(num_temp));
	AllocateSpace(num_temp);
	if (num_temp>0)
		file.read((char*)pdata, num_temp*sizeof(char));
}

////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
MCMC_PixelResp::MCMC_PixelResp()
{
	presps = NULL;
	num_value = 0;
}

MCMC_PixelResp::~MCMC_PixelResp()
{
	Free();
}

void MCMC_PixelResp::Free(void)
{
	if (presps != NULL)
		delete []presps;
	presps = NULL;
	num_value = 0;
}

void MCMC_PixelResp::AllocateSpace(const int size)
{
	if (size != num_value)
	{
		Free();

		num_value = size;
		if (num_value>0)
			presps = new MCMC_RespChange[size];
	}
}

void MCMC_PixelResp::SaveToFile(ofstream &file)
{
	file.write((char*)&num_value, sizeof(num_value));
	for (int i=0; i<num_value; i++)
		presps[i].SaveToFile(file);
}

void MCMC_PixelResp::ReadFromFile(ifstream &file)
{
	int num_temp;

	file.read((char*)&num_temp, sizeof(int));
	AllocateSpace(num_temp);

	for (int i=0; i<num_value; i++)
		presps[i].ReadFromFile(file);
}

////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
MCMC_PLE_Resp::MCMC_PLE_Resp()
{
	p_filter_bank = NULL;
	bRememberAsHist = false;
}

MCMC_PLE_Resp::~MCMC_PLE_Resp()
{
	p_filter_bank = NULL;
}

void MCMC_PLE_Resp::SetFiter(MCMC_FilterBank &filter_bank, const int ifilter,
							 const int rows, const int cols,
							 const double low_input, const double high_input, const double delta_input,
							 const double low_output, const double high_output, const double delta_output)
{
	p_filter_bank = &filter_bank;
	nfilter = ifilter;

	mx_filter_resp.SetDimension(rows,cols);

	bRememberAsHist = false;

	// allocate space
	hist_ref_input.Set(low_input, high_input, delta_input);
	hist_ref_output.Set(low_output, high_output, delta_output);

	MCMC_PixelResp	*p_pixel;
	MCMC_RespChange	*p_resp;

	for (int j=0; j<mx_filter_resp.rows(); j++)
		for (int i=0; i<mx_filter_resp.cols(); i++)
		{
			p_pixel = &(mx_filter_resp(j,i));
			p_pixel->AllocateSpace(hist_ref_input.count);

			for (int k=0; k<hist_ref_input.count; k++)
			{
				p_resp = &(p_pixel->presps[k]);
				if (bRememberAsHist)
					p_resp->AllocateSpace(hist_ref_output.count+2);	// save as histogram
				else
					p_resp->AllocateSpace(p_filter_bank->p_filters[nfilter].NumSparse());		// save as resposes
			}
		}
}

void MCMC_PLE_Resp::ComputeResponse(const int row_top, const int row_bottom, const int col_left, const int col_right)
{
	MCMC_PixelResp	*p_pixel;
	MCMC_RespChange	*p_resp;
	int				m,n,j,i,k,t;
	double			value;
	mcmcVector<double>	vc_resp;

	// bouded rectangle
	rect_bound.left		= col_left;
	rect_bound.right	= col_right;
	rect_bound.top		= row_top;
	rect_bound.bottom	= row_bottom;

	for (j=rect_bound.top; j<=rect_bound.bottom; j++)
		for (i=rect_bound.left; i<=rect_bound.right; i++)
		{
			m = j-row_top;
			n = i-col_left;

			p_pixel = &(mx_filter_resp(m,n));

			for (k=0; k<hist_ref_input.count; k++)
			{
				value = hist_ref_input.GetValueIndex(k);

				p_resp = &(p_pixel->presps[k]);
				if (bRememberAsHist)
				{
					assert(false);
/*
					hist_ref_output.Init(0.0);
					p_filter_bank->NewFilterRespAtOneSite(hist_ref_output,nfilter,j,i,value);

					for (t=0; t<hist_ref_output.count; t++)
						p_resp->pdata[t] = hist_ref_output.pbins[t];
					// values that are out of bounds
					p_resp->pdata[hist_ref_output.count] = hist_ref_output.outof_lower_bound;
					p_resp->pdata[hist_ref_output.count+1] = hist_ref_output.outof_upper_bound;
*/
				}
				else
				{
					p_filter_bank->NewFilterRespAtOneSite(vc_resp,nfilter,j,i,value);
					for (t=0; t<p_resp->num_resp; t++)
						p_resp->pdata[t] = hist_ref_output.GetIndex(vc_resp(t));
				}
			}
		}
}

void MCMC_PLE_Resp::SaveToFile(ofstream &file)
{
	file.write((char*)&rect_bound.top, sizeof(rect_bound.top));
	file.write((char*)&rect_bound.bottom, sizeof(rect_bound.bottom));
	file.write((char*)&rect_bound.left, sizeof(rect_bound.left));
	file.write((char*)&rect_bound.right, sizeof(rect_bound.right));

	for (int j=0; j<rect_bound.Height(); j++)
		for (int i=0; i<rect_bound.Width(); i++)
			mx_filter_resp(j,i).SaveToFile(file);
}

void MCMC_PLE_Resp::ReadFromFile(ifstream &file)
{
	file.read((char*)&rect_bound.top, sizeof(rect_bound.top));
	file.read((char*)&rect_bound.bottom, sizeof(rect_bound.bottom));
	file.read((char*)&rect_bound.left, sizeof(rect_bound.left));
	file.read((char*)&rect_bound.right, sizeof(rect_bound.right));

	for (int j=0; j<rect_bound.Height(); j++)
		for (int i=0; i<rect_bound.Width(); i++)
			mx_filter_resp(j,i).ReadFromFile(file);
}

////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////
FRAME_One_Channel::FRAME_One_Channel()
{
	p_ple_resp = NULL;

	rf_rand_beta.Set(-1.0,1.0,0.1);
}

FRAME_One_Channel::~FRAME_One_Channel()
{
	if (p_ple_resp != NULL)
		delete []p_ple_resp;
	p_ple_resp = NULL;
}

void FRAME_One_Channel::Set(MCMC_FilterBank &filter_bank, const bool bComputeInAdvance,
							const int biggest_row_size, const int biggest_col_size)
{
	int i;

	// initialization
	pfilter_bank = &filter_bank;
	num_filter = filter_bank.NumFilter();
	vc_histogram_obs.SetDimension(num_filter,1);
	vc_histogram_syn.SetDimension(num_filter,1);
	vc_histogram_syn_temp.SetDimension(num_filter,1);
	vc_B.SetDimension(num_filter,1);

	for (i=0; i<NumFilter(); i++)
	{
		vc_histogram_obs(i) = filter_bank.p_histogram_obs[i];
		vc_histogram_syn(i).Set(filter_bank.p_min_value[i],filter_bank.p_max_value[i],filter_bank.p_deltas[i]);
		vc_histogram_syn(i).Init(0.0);
		vc_histogram_syn_temp(i).Set(filter_bank.p_min_value[i],filter_bank.p_max_value[i],filter_bank.p_deltas[i]);
		vc_histogram_syn_temp(i).Init(0.0);

		vc_B(i).Set(filter_bank.p_min_value[i],filter_bank.p_max_value[i],filter_bank.p_deltas[i]);
		vc_B(i).Init(1.0);

		// parameters for the values that are out of bounds
		vc_B(i).SetB4OutLowerBound(90.0);
		vc_B(i).SetB4OutUpperBound(90.0);
	}
	if (NumFilter()==1)
	{
		// if only delta filter is used then use the histogram directly
		for (i=0; i<vc_B(0).count; i++)
		{
			if (vc_histogram_obs(0).pbins[i]<1e-10)
				vc_B(0).pbins[i] = 20;
			else
				vc_B(0).pbins[i] = -log(vc_histogram_obs(0).pbins[i]);
		}
	}
	else
	{
		for (int j=0; j<NumFilter(); j++)
		{
			bool bset_delta=false;
			if (j==0)
			{
				for (i=0; i<vc_B(j).count; i++)
				{
					if (vc_histogram_obs(j).pbins[i]>0.70)
						bset_delta = true;
				}
				bset_delta = true;
			}

			for (i=0; i<vc_B(j).count; i++)
			{
				if (vc_histogram_obs(j).pbins[i]<1e-10)
					vc_B(j).pbins[i] = 20;
				else
				{
					if (!bset_delta)
						vc_B(j).pbins[i] = 1.0;
					else
						vc_B(j).pbins[i] = -log(vc_histogram_obs(0).pbins[i]);
				}
			}
		}
	}

	low_value	= filter_bank.p_min_value[0];
	high_value	= filter_bank.p_max_value[0];
	delta		= filter_bank.p_deltas[0];

	histogram_ref_input.Set(low_value,high_value,delta);

	// saving energies temporarily
	vc_energy.SetRowVector(NumValue());
	vc_energy.InitValue(1000.0);

	if (bComputeInAdvance)
	{
		if (p_ple_resp != NULL)
			delete []p_ple_resp;
		p_ple_resp = new MCMC_PLE_Resp[NumFilter()];
		for (i=0; i<NumFilter(); i++)
			p_ple_resp[i].SetFiter(filter_bank, i, biggest_row_size, biggest_col_size,
								   GetLow(), GetHigh(), GetDelta(),
								   filter_bank.p_min_value[i], filter_bank.p_max_value[i], filter_bank.p_deltas[i]);
	}
}

void FRAME_One_Channel::ComputeResponse(const int row_top, const int row_bottom, const int col_left, const int col_right)
{
	for (int i=0; i<NumFilter(); i++)
		p_ple_resp[i].ComputeResponse(row_top,row_bottom,col_left,col_right);
}

void FRAME_One_Channel::SaveRespHistToFile(ofstream &file)
{
	for (int i=0; i<num_filter; i++)
		p_ple_resp[i].SaveToFile(file);
}

void FRAME_One_Channel::ReadRespHistFromFile(ifstream &file)
{
	for (int i=0; i<num_filter; i++)
		p_ple_resp[i].ReadFromFile(file);
}

void FRAME_One_Channel::FilteringOnNewImage(Raster<double> &image)
{
	if (pfilter_bank != NULL)
	{
		pfilter_bank->ras_image = image;
		pfilter_bank->DoFiltering();
	}
}

bool FRAME_One_Channel::UpdateBeta(const double d)
{
	bool bConverge=true;
	mcmcVector<double> vc_dist;
	int				i,j;

	for (i=0; i<num_filter; i++)
	{
		vc_histogram_syn(i).Normalize();
		vc_dist.SetRowVector(vc_histogram_obs(i).count);
		for (j=0; j<vc_dist.size(); j++)
			vc_dist(j) = vc_histogram_obs(i).pbins[j]-vc_histogram_syn(i).pbins[j];

		if (bConverge && vc_dist.Norm(1.0)>=0.07)
			bConverge = false;
		for (j=0; j<vc_dist.size(); j++)
			vc_B(i).pbins[j] = vc_B(i).pbins[j]-vc_dist(j)*d;
	}
	return bConverge;
}

double FRAME_One_Channel::EnergyFromHist(const int row, const int col,const int index_value)
{
	double			value;
	double			energy=0.0;
	int				i,j;
	MCMC_PixelResp	*p_pixel_resp;
	MCMC_RespChange	*p_resp;
	Gibbs_B			*p_B;

	value = histogram_ref_input.GetValueIndex(index_value);

	for (i=0; i<num_filter; i++)
	{
		p_B = &(vc_B(i));
		// get pointer to the responses at pixel (row,col)
		p_pixel_resp = &(p_ple_resp[i].mx_filter_resp(row,col));
		// get all the responses w.r.t. 'value' at pixel (row,col)
		p_resp = &(p_pixel_resp->presps[index_value]);
		if (p_ple_resp[i].AsHistogram())
		{
assert(false);
			for (j=0; j<p_B->count; j++)
				energy += p_B->pbins[j]*p_resp->pdata[j];
			// out of bounds
			energy += p_B->outof_lower_bound*p_resp->pdata[p_B->count];
			energy += p_B->outof_upper_bound*p_resp->pdata[p_B->count+1];
		}
		else
		{
			for (j=0; j<p_resp->num_resp; j++)
			{
				if (p_resp->pdata[j]<0)
					energy += p_B->outof_lower_bound;
				else if(p_resp->pdata[j]>=p_B->count)
					energy += p_B->outof_upper_bound;
				else
					energy += p_B->pbins[p_resp->pdata[j]];
			}
		}
	}
	return energy;
}

void FRAME_One_Channel::UpdateHistogramSyn(void)
{
	int		i,j;
	int	size;

	for (i=0; i<num_filter; i++)
	{
		vc_histogram_syn_temp(i).Normalize();
		size = vc_histogram_syn_temp(i).count;
		for (j=0; j<size; j++)
			vc_histogram_syn(i).pbins[j] += vc_histogram_syn_temp(i).pbins[j];
		vc_histogram_syn(i).outof_lower_bound += vc_histogram_syn_temp(i).outof_lower_bound;
		vc_histogram_syn(i).outof_upper_bound += vc_histogram_syn_temp(i).outof_upper_bound;
	}
}

void FRAME_One_Channel::UpdateHistogramSynTemp(const double prob, const double value, const int row, const int col)
{
	int				i,j,index_value;
	MCMC_PixelResp	*p_pixel_resp;
	MCMC_RespChange	*p_resp;


	if (prob > 0.0)
	{
		index_value = GetIndex(value);

		for (i=0; i<num_filter; i++)
		{
			// get pointer to the responses at pixel (row,col)
			p_pixel_resp = &(p_ple_resp[i].mx_filter_resp(row,col));
			// get all the responses w.r.t. 'value' at pixel (row,col)
			p_resp = &(p_pixel_resp->presps[index_value]);

			Histogram	*phist=&(vc_histogram_syn_temp(i));

			if (p_ple_resp[i].AsHistogram())
			{
assert(false);
				for (j=0; j<phist->count; j++)
					phist->pbins[j] += p_resp->pdata[j]*prob;
				phist->outof_lower_bound += p_resp->pdata[phist->count]*prob;
				phist->outof_upper_bound += p_resp->pdata[phist->count+1]*prob;
			}
			else
			{
				for (j=0; j<p_resp->num_resp; j++)
				{
					if (p_resp->pdata[j]>=0 && p_resp->pdata[j]<phist->count)
						phist->pbins[p_resp->pdata[j]] += prob;
				}
			}
		}
	}
}

void FRAME_One_Channel::ClearHistogramSyn(void)
{
	for (int i=0; i<num_filter; i++)
		vc_histogram_syn(i).Init(0.0);
}

void FRAME_One_Channel::ClearHistogramSynTemp(void)
{
	for (int i=0; i<num_filter; i++)
		vc_histogram_syn_temp(i).Init(0.0);
}

void FRAME_One_Channel::ComputeEnergyAtEachValue(const int j, const int i)
{
	double			value;
	int				k,size,t,nfilter;
	Vector_Double	vc_resp;
	Gibbs_B			*p_B;

	for (k=0; k<NumValue(); k++)
	{
		value = histogram_ref_input.GetValueIndex(k);

		vc_energy(k) = 0.0;

		for (nfilter=0; nfilter<NumFilter(); nfilter++)
		{
			pfilter_bank->NewFilterRespAtOneSite(vc_resp,nfilter,j,i,value);
			p_B = &(vc_B(nfilter));

			size = vc_resp.size();
			for (t=0; t<size; t++)
				vc_energy(k) += p_B->GetBin(vc_resp(t));
		}
	}
}

// from the pre-computed histogram
void FRAME_One_Channel::SaveTempEnergy(const int row, const int col)
{
	for (int k=0; k<NumValue(); k++)
		vc_energy(k) = EnergyFromHist(row,col,k);
}
